// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_BANDMATRIX_H
#define EIGEN_BANDMATRIX_H

namespace Eigen {

namespace internal {

template<typename Derived>
class BandMatrixBase : public EigenBase<Derived>
{
  public:
	enum
	{
		Flags = internal::traits<Derived>::Flags,
		CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
		RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
		ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
		MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
		MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
		Supers = internal::traits<Derived>::Supers,
		Subs = internal::traits<Derived>::Subs,
		Options = internal::traits<Derived>::Options
	};
	typedef typename internal::traits<Derived>::Scalar Scalar;
	typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime> DenseMatrixType;
	typedef typename DenseMatrixType::StorageIndex StorageIndex;
	typedef typename internal::traits<Derived>::CoefficientsType CoefficientsType;
	typedef EigenBase<Derived> Base;

  protected:
	enum
	{
		DataRowsAtCompileTime = ((Supers != Dynamic) && (Subs != Dynamic)) ? 1 + Supers + Subs : Dynamic,
		SizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime)
	};

  public:
	using Base::cols;
	using Base::derived;
	using Base::rows;

	/** \returns the number of super diagonals */
	inline Index supers() const { return derived().supers(); }

	/** \returns the number of sub diagonals */
	inline Index subs() const { return derived().subs(); }

	/** \returns an expression of the underlying coefficient matrix */
	inline const CoefficientsType& coeffs() const { return derived().coeffs(); }

	/** \returns an expression of the underlying coefficient matrix */
	inline CoefficientsType& coeffs() { return derived().coeffs(); }

	/** \returns a vector expression of the \a i -th column,
	 * only the meaningful part is returned.
	 * \warning the internal storage must be column major. */
	inline Block<CoefficientsType, Dynamic, 1> col(Index i)
	{
		EIGEN_STATIC_ASSERT((int(Options) & int(RowMajor)) == 0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
		Index start = 0;
		Index len = coeffs().rows();
		if (i <= supers()) {
			start = supers() - i;
			len = (std::min)(rows(), std::max<Index>(0, coeffs().rows() - (supers() - i)));
		} else if (i >= rows() - subs())
			len = std::max<Index>(0, coeffs().rows() - (i + 1 - rows() + subs()));
		return Block<CoefficientsType, Dynamic, 1>(coeffs(), start, i, len, 1);
	}

	/** \returns a vector expression of the main diagonal */
	inline Block<CoefficientsType, 1, SizeAtCompileTime> diagonal()
	{
		return Block<CoefficientsType, 1, SizeAtCompileTime>(coeffs(), supers(), 0, 1, (std::min)(rows(), cols()));
	}

	/** \returns a vector expression of the main diagonal (const version) */
	inline const Block<const CoefficientsType, 1, SizeAtCompileTime> diagonal() const
	{
		return Block<const CoefficientsType, 1, SizeAtCompileTime>(
			coeffs(), supers(), 0, 1, (std::min)(rows(), cols()));
	}

	template<int Index>
	struct DiagonalIntReturnType
	{
		enum
		{
			ReturnOpposite =
				(int(Options) & int(SelfAdjoint)) && (((Index) > 0 && Supers == 0) || ((Index) < 0 && Subs == 0)),
			Conjugate = ReturnOpposite && NumTraits<Scalar>::IsComplex,
			ActualIndex = ReturnOpposite ? -Index : Index,
			DiagonalSize =
				(RowsAtCompileTime == Dynamic || ColsAtCompileTime == Dynamic)
					? Dynamic
					: (ActualIndex < 0
						   ? EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime, RowsAtCompileTime + ActualIndex)
						   : EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime - ActualIndex))
		};
		typedef Block<CoefficientsType, 1, DiagonalSize> BuildType;
		typedef typename internal::conditional<Conjugate,
											   CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, BuildType>,
											   BuildType>::type Type;
	};

	/** \returns a vector expression of the \a N -th sub or super diagonal */
	template<int N>
	inline typename DiagonalIntReturnType<N>::Type diagonal()
	{
		return typename DiagonalIntReturnType<N>::BuildType(
			coeffs(), supers() - N, (std::max)(0, N), 1, diagonalLength(N));
	}

	/** \returns a vector expression of the \a N -th sub or super diagonal */
	template<int N>
	inline const typename DiagonalIntReturnType<N>::Type diagonal() const
	{
		return typename DiagonalIntReturnType<N>::BuildType(
			coeffs(), supers() - N, (std::max)(0, N), 1, diagonalLength(N));
	}

	/** \returns a vector expression of the \a i -th sub or super diagonal */
	inline Block<CoefficientsType, 1, Dynamic> diagonal(Index i)
	{
		eigen_assert((i < 0 && -i <= subs()) || (i >= 0 && i <= supers()));
		return Block<CoefficientsType, 1, Dynamic>(coeffs(), supers() - i, std::max<Index>(0, i), 1, diagonalLength(i));
	}

	/** \returns a vector expression of the \a i -th sub or super diagonal */
	inline const Block<const CoefficientsType, 1, Dynamic> diagonal(Index i) const
	{
		eigen_assert((i < 0 && -i <= subs()) || (i >= 0 && i <= supers()));
		return Block<const CoefficientsType, 1, Dynamic>(
			coeffs(), supers() - i, std::max<Index>(0, i), 1, diagonalLength(i));
	}

	template<typename Dest>
	inline void evalTo(Dest& dst) const
	{
		dst.resize(rows(), cols());
		dst.setZero();
		dst.diagonal() = diagonal();
		for (Index i = 1; i <= supers(); ++i)
			dst.diagonal(i) = diagonal(i);
		for (Index i = 1; i <= subs(); ++i)
			dst.diagonal(-i) = diagonal(-i);
	}

	DenseMatrixType toDenseMatrix() const
	{
		DenseMatrixType res(rows(), cols());
		evalTo(res);
		return res;
	}

  protected:
	inline Index diagonalLength(Index i) const
	{
		return i < 0 ? (std::min)(cols(), rows() + i) : (std::min)(rows(), cols() - i);
	}
};

/**
 * \class BandMatrix
 * \ingroup Core_Module
 *
 * \brief Represents a rectangular matrix with a banded storage
 *
 * \tparam _Scalar Numeric type, i.e. float, double, int
 * \tparam _Rows Number of rows, or \b Dynamic
 * \tparam _Cols Number of columns, or \b Dynamic
 * \tparam _Supers Number of super diagonal
 * \tparam _Subs Number of sub diagonal
 * \tparam _Options A combination of either \b #RowMajor or \b #ColMajor, and of \b #SelfAdjoint
 *                  The former controls \ref TopicStorageOrders "storage order", and defaults to
 *                  column-major. The latter controls whether the matrix represents a selfadjoint
 *                  matrix in which case either Supers of Subs have to be null.
 *
 * \sa class TridiagonalMatrix
 */

template<typename _Scalar, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
struct traits<BandMatrix<_Scalar, _Rows, _Cols, _Supers, _Subs, _Options>>
{
	typedef _Scalar Scalar;
	typedef Dense StorageKind;
	typedef Eigen::Index StorageIndex;
	enum
	{
		CoeffReadCost = NumTraits<Scalar>::ReadCost,
		RowsAtCompileTime = _Rows,
		ColsAtCompileTime = _Cols,
		MaxRowsAtCompileTime = _Rows,
		MaxColsAtCompileTime = _Cols,
		Flags = LvalueBit,
		Supers = _Supers,
		Subs = _Subs,
		Options = _Options,
		DataRowsAtCompileTime = ((Supers != Dynamic) && (Subs != Dynamic)) ? 1 + Supers + Subs : Dynamic
	};
	typedef Matrix<Scalar, DataRowsAtCompileTime, ColsAtCompileTime, int(Options) & int(RowMajor) ? RowMajor : ColMajor>
		CoefficientsType;
};

template<typename _Scalar, int Rows, int Cols, int Supers, int Subs, int Options>
class BandMatrix : public BandMatrixBase<BandMatrix<_Scalar, Rows, Cols, Supers, Subs, Options>>
{
  public:
	typedef typename internal::traits<BandMatrix>::Scalar Scalar;
	typedef typename internal::traits<BandMatrix>::StorageIndex StorageIndex;
	typedef typename internal::traits<BandMatrix>::CoefficientsType CoefficientsType;

	explicit inline BandMatrix(Index rows = Rows, Index cols = Cols, Index supers = Supers, Index subs = Subs)
		: m_coeffs(1 + supers + subs, cols)
		, m_rows(rows)
		, m_supers(supers)
		, m_subs(subs)
	{
	}

	/** \returns the number of columns */
	inline EIGEN_CONSTEXPR Index rows() const { return m_rows.value(); }

	/** \returns the number of rows */
	inline EIGEN_CONSTEXPR Index cols() const { return m_coeffs.cols(); }

	/** \returns the number of super diagonals */
	inline EIGEN_CONSTEXPR Index supers() const { return m_supers.value(); }

	/** \returns the number of sub diagonals */
	inline EIGEN_CONSTEXPR Index subs() const { return m_subs.value(); }

	inline const CoefficientsType& coeffs() const { return m_coeffs; }
	inline CoefficientsType& coeffs() { return m_coeffs; }

  protected:
	CoefficientsType m_coeffs;
	internal::variable_if_dynamic<Index, Rows> m_rows;
	internal::variable_if_dynamic<Index, Supers> m_supers;
	internal::variable_if_dynamic<Index, Subs> m_subs;
};

template<typename _CoefficientsType, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
class BandMatrixWrapper;

template<typename _CoefficientsType, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
struct traits<BandMatrixWrapper<_CoefficientsType, _Rows, _Cols, _Supers, _Subs, _Options>>
{
	typedef typename _CoefficientsType::Scalar Scalar;
	typedef typename _CoefficientsType::StorageKind StorageKind;
	typedef typename _CoefficientsType::StorageIndex StorageIndex;
	enum
	{
		CoeffReadCost = internal::traits<_CoefficientsType>::CoeffReadCost,
		RowsAtCompileTime = _Rows,
		ColsAtCompileTime = _Cols,
		MaxRowsAtCompileTime = _Rows,
		MaxColsAtCompileTime = _Cols,
		Flags = LvalueBit,
		Supers = _Supers,
		Subs = _Subs,
		Options = _Options,
		DataRowsAtCompileTime = ((Supers != Dynamic) && (Subs != Dynamic)) ? 1 + Supers + Subs : Dynamic
	};
	typedef _CoefficientsType CoefficientsType;
};

template<typename _CoefficientsType, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
class BandMatrixWrapper
	: public BandMatrixBase<BandMatrixWrapper<_CoefficientsType, _Rows, _Cols, _Supers, _Subs, _Options>>
{
  public:
	typedef typename internal::traits<BandMatrixWrapper>::Scalar Scalar;
	typedef typename internal::traits<BandMatrixWrapper>::CoefficientsType CoefficientsType;
	typedef typename internal::traits<BandMatrixWrapper>::StorageIndex StorageIndex;

	explicit inline BandMatrixWrapper(const CoefficientsType& coeffs,
									  Index rows = _Rows,
									  Index cols = _Cols,
									  Index supers = _Supers,
									  Index subs = _Subs)
		: m_coeffs(coeffs)
		, m_rows(rows)
		, m_supers(supers)
		, m_subs(subs)
	{
		EIGEN_UNUSED_VARIABLE(cols);
		// internal::assert(coeffs.cols()==cols() && (supers()+subs()+1)==coeffs.rows());
	}

	/** \returns the number of columns */
	inline EIGEN_CONSTEXPR Index rows() const { return m_rows.value(); }

	/** \returns the number of rows */
	inline EIGEN_CONSTEXPR Index cols() const { return m_coeffs.cols(); }

	/** \returns the number of super diagonals */
	inline EIGEN_CONSTEXPR Index supers() const { return m_supers.value(); }

	/** \returns the number of sub diagonals */
	inline EIGEN_CONSTEXPR Index subs() const { return m_subs.value(); }

	inline const CoefficientsType& coeffs() const { return m_coeffs; }

  protected:
	const CoefficientsType& m_coeffs;
	internal::variable_if_dynamic<Index, _Rows> m_rows;
	internal::variable_if_dynamic<Index, _Supers> m_supers;
	internal::variable_if_dynamic<Index, _Subs> m_subs;
};

/**
 * \class TridiagonalMatrix
 * \ingroup Core_Module
 *
 * \brief Represents a tridiagonal matrix with a compact banded storage
 *
 * \tparam Scalar Numeric type, i.e. float, double, int
 * \tparam Size Number of rows and cols, or \b Dynamic
 * \tparam Options Can be 0 or \b SelfAdjoint
 *
 * \sa class BandMatrix
 */
template<typename Scalar, int Size, int Options>
class TridiagonalMatrix : public BandMatrix<Scalar, Size, Size, Options & SelfAdjoint ? 0 : 1, 1, Options | RowMajor>
{
	typedef BandMatrix<Scalar, Size, Size, Options & SelfAdjoint ? 0 : 1, 1, Options | RowMajor> Base;
	typedef typename Base::StorageIndex StorageIndex;

  public:
	explicit TridiagonalMatrix(Index size = Size)
		: Base(size, size, Options & SelfAdjoint ? 0 : 1, 1)
	{
	}

	inline typename Base::template DiagonalIntReturnType<1>::Type super() { return Base::template diagonal<1>(); }
	inline const typename Base::template DiagonalIntReturnType<1>::Type super() const
	{
		return Base::template diagonal<1>();
	}
	inline typename Base::template DiagonalIntReturnType<-1>::Type sub() { return Base::template diagonal<-1>(); }
	inline const typename Base::template DiagonalIntReturnType<-1>::Type sub() const
	{
		return Base::template diagonal<-1>();
	}

  protected:
};

struct BandShape
{};

template<typename _Scalar, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
struct evaluator_traits<BandMatrix<_Scalar, _Rows, _Cols, _Supers, _Subs, _Options>>
	: public evaluator_traits_base<BandMatrix<_Scalar, _Rows, _Cols, _Supers, _Subs, _Options>>
{
	typedef BandShape Shape;
};

template<typename _CoefficientsType, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
struct evaluator_traits<BandMatrixWrapper<_CoefficientsType, _Rows, _Cols, _Supers, _Subs, _Options>>
	: public evaluator_traits_base<BandMatrixWrapper<_CoefficientsType, _Rows, _Cols, _Supers, _Subs, _Options>>
{
	typedef BandShape Shape;
};

template<>
struct AssignmentKind<DenseShape, BandShape>
{
	typedef EigenBase2EigenBase Kind;
};

} // end namespace internal

} // end namespace Eigen

#endif // EIGEN_BANDMATRIX_H
